Analysis

BUSINESS QUESTIONS: Note: ignore the fact there is no question 1, this is due to a formatting error 2. a) Is there are a trend in injuries By Region, and is there different peaks of times of year per region? b) Analysis on the types of animals that are injured, this also by Region – is there a species that is more liable to injury in certain regions? - (cat/dog by region) c) What is the outcome? Does this differ by region?

1. Total call volume for complaint calls: How has this trended over time?
2. Is there a particular animal being called about the most?
3. Do particular suburbs have different type of complaint calls? Do they call about different animals? ((MAKE A LEAFLET MAP FOR THIS!))

4. Business Intelligence – using the insights you have found, can you predict how this might look for the upcoming year?

source("cleaning_script.R")

Rows: 31330 Columns: 7── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr (6): nature, animal_type, category, suburb, date_range, city
lgl (1): responsible_office
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Warning: Expected 2 pieces. Missing pieces filled with `NA` in 5643 rows [7, 11, 12, 18, 19, 25, 27, 29, 30, 31, 39, 44, 52, 54, 66, 67, 74, 75, 78, 81, ...].Rows: 42413 Columns: 5── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): Animal Type, Complaint Type, Date Received, Suburb, Electoral Division
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 664 Columns: 12── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): animal_type, outcome
dbl (10): year, ACT, NSW, NT, QLD, SA, TAS, VIC, WA, Total
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Note; summer in Australia seasons: Summer: December - February Autumn: March - May Winter: June - August Spring: September - November

also: The ‘wet season’ in Australia’s North: November - April

https://www.mdpi.com/2076-2615/8/7/100 will useful reading for later, talk about the need to reduce euthanasia or similar and the effect this has on people

Intro: Introduce the point of the talk, talk about Australian RSPCA, talk about how the data was gathered (using the information from the websites), the PURPOSE of this investigation (which is to help the RSPCA know which areas/animals to focus their efforts on), as I’m introducing the datasets I can introduce the two different cities. First, we can talk a short bit about Australia as a whole, it’s climate, the kind of animals etc. The point here is to really set the scene before diving too deep into facts and figures, as especially a non-technical audience this will help keep them engaged and make a more holistic presentation. In my opinion, it’s always good to zoom out and see the big picture, rather than getting lost in the myopia of some csv files.

Townsville Intro: Townsville is a city on the north-eastern coast of Queensland, Australia. With a population of 180,820 as of June 2018, it is the largest settlement in North Queensland; it is unofficially considered its capital. [note: put a map of Australia with Queensland and Townsville highlighted here. Talk a little bit about the population density, urbanisation, climate, types of animals that are common here etc. Show some photos of the area too]

Brisbane Intro: Exact same as above

3.

1. Total call volume for complaint calls: How has this trended over time?

First, let’s look at the Townsville animal complaints.

animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date_received, y = count)) +
  geom_line() +
  scale_x_date(date_breaks = "6 months", date_labels = "%b-%y")

This shows some seasonality and also an increase then a decline. Each summer (in December) the calls are much lower, rising again each Winter. Using geom_smooth, we get:

animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date_received, y = count)) +
  geom_smooth() +
  geom_point() +
  scale_x_date(date_breaks = "1 year", date_labels = "%y")

This shows a bit more clearly the general trend of call volume. From 2014 it steadily rises, peaking in 2017. Afterwards, it steadily declines to only slightly higher than where it started. It would be very difficult to say whether this trend will continue downward, go upward or stay relatively flat.

Note; if we have time or if it’s helpful, we will fix this graph so that quarters are properly displayed.

brisbane_complaints %>% 
  group_by(date) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date, y = count)) +
  geom_point() +
  geom_line() +
  scale_x_date(date_breaks = "3 months", date_labels = "%y")

Again, we see the seasonality of winter having more calls. The fact this is in both Brisbane and Townsville suggests a fairly general trend.

brisbane_complaints %>% 
  group_by(date) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date, y = count)) +
  geom_point() +
  geom_smooth() +
  scale_x_date(date_breaks = "1 year", date_labels = "20%y")

Other than declining a little over 2016, the number of calls sees a slow but steady increase towards 2020, being thousands more than it was in 2016 and 2017. So, the general trend is that the RSPCA are getting more complaints as time goes on. Now, this does not necessarily mean the line will continue to go up. We are also missing Q3 from 2016 so this skews the curve a little

3.

2. Is there a particular animal being called about the most?

animal_complaints %>% 
  group_by(animal_type) %>% 
  ggplot(aes(x = animal_type)) +
  geom_bar()

The number of calls about dogs dwarf those about cats hugely. Let’s look at specific numbers:

animal_complaints %>% 
  group_by(animal_type) %>%
  count() %>% 
  summarise(n, 4094 / 38319)

So cats account for only 10% of the calls dogs do for Townsville! Let’s look at Brisbane:

brisbane_complaints %>% 
  group_by(type_of_animal) %>%
  count()

brisbane_complaints %>% 
  group_by(type_of_animal) %>%
  ggplot(aes(x = type_of_animal)) +
  geom_bar()

We have a lot more animal types here that we don’t in the previous data. We also have many calls about Attacks with no animal specified, many of which it is very likely they were involving dogs. However, we can’t say this for sure.

So that makes

4745 / 13334  # cats are 35% of the calls compared to dogs

[1] 0.3558572

# all other animals divided by dogs, leaving out unspecified (which we believe potentially contain a high proportion of dogs)
9457 / 13349

[1] 0.7084426

Most of the other animals have very small counts, interestingly foxes seem to make up a decent proportion of the calls regarding wild animals.

To answer the question however, it’s mainly dogs and cats, and especially dogs. Dogs are being called about more than any of the other animals combined (if we leave Unspecified to the side)

(Attack refers to the initial description of the complaint)

3.

3. Do particular suburbs have different type of complaint calls? Do they call about different animals?

OK, so this one is difficult simply for the fact there is a huge amount of suburbs.

brisbane_complaints %>% 
  group_by(suburb) %>% 
  count()

# There is 192 suburbs! Certainly using a fill on a graph is not gonna work, neither is putting them on the x-axis on a bar graph.

animal_complaints %>% 
  group_by(suburb) %>% 
  count()

# 85 suburbs for Townsville

animal_complaints %>% 
  group_by(electoral_division) %>% 
  count()

# also 11 electoral divisions. Could we look at suburbs one electoral division at a time? Possibly

Idea: What about using leaflet to visualise types of complaint calls on a map?

This is more feasible with the Townsville datasets, as it has less suburbs and only 6 types of complaints. It may be necessary to try and wrangle the Brisbane data a little further to narrow down categories (for both type_of_animal and complaint_type)

FOR TOWNSVILLE WE HAVE: 85 different suburbs 6 complaint types 2 animal types

We want to break it down by suburb and complaint type, and then by suburb and animal type

animal_complaints %>% 
  ggplot(aes(x = suburb, fill = animal_type)) +
  geom_bar(position = "fill") +
  coord_flip() +
  scale_x_discrete(guide = guide_axis(n.dodge = 5))

animal_complaints %>% 
  ggplot(aes(x = suburb, fill =complaint_type)) +
  geom_bar(position = "fill") 

  #coord_flip() +
  #scale_x_discrete(guide = guide_axis(n.dodge = 5))

The challenge again, is having too many suburbs that we can’t get any useful information out of the graphs

animal_complaints %>% 
  group_by(suburb) %>% 
  count(sort = TRUE)

# let's drop any suburbs with less than 200 cases

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), animal_type) %>% 
  filter(count >= 500) %>% 
  ggplot(aes(x = suburb, y = count, fill = animal_type)) +
  geom_col() +
  coord_flip()

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

The highest count by quite a large margin is unallocated to a specific suburb.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), animal_type) %>% 
  filter(count >= 500) %>% 
  ggplot(aes(x = suburb, y = count, fill = animal_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

At a glance, there’s no big difference on animal types. All the suburbs have a large majority of dogs. Could we do a hypothesis test to see if there is a statistically significant difference?

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), animal_type) %>% 
  filter(count < 500 & count > 100) %>% 
  ggplot(aes(x = suburb, y = count, fill = animal_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

Looking at the lower end of the population, we see an outlier. Townsville City has way more cats than any of the other suburbs. In fact, it’s almost 50/50!

Now still looking at the Townsville dataset, we’ll break it down by complaint type:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500 & count <= 4000) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

Now for the lower end of the count:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 50) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

I’m not sure much can be gleamed with so many complaint_types. Let’s try focusing more specifically.

One thing to note here, is the vast variation in tolerance for noise. Most categories are consistent except this one. Cluden for example has a small percentage of noise complaints,while Bohle Plains has a huge percentage. This requires more investigation to figure out the root cause of this:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Noise" | (complaint_type == "Attack")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Noise" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Wandering" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Enclosure" | (complaint_type == "Noise")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

For the less than 500 greater than 50 groups:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>% 
  filter(complaint_type == "Noise" | (complaint_type == "Attack")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>%  
  filter(complaint_type == "Noise" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>% 
  filter(complaint_type == "Wandering" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>% 
  #filter(complaint_type == "Enclosure" | (complaint_type == "Noise")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() +
  facet_wrap(~ suburb)

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type, animal_type) %>% 
  filter(count < 500 & count > 100) %>% 
  ggplot(aes(x = complaint_type, y = count, fill = complaint_type)) +
  geom_col() +
  facet_wrap(~ suburb)

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

From this graph, we can see that Hyde Park has a disproportionate amount of private impounds, while Bohle plains has more noise complaints.

Potential hypothesis test:

That the level of private impounds in Hyde Park being greater is statistically significant.

That the level of noise in Bohle Plains being greater is statistically significant.

animal_complaints %>% 
  group_by(suburb) %>% 
  count() %>% 
  filter(n < 500 & n > 100)

Let’s look at the suburbs with much higher complaints now:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type, animal_type) %>% 
  filter(count >= 500 & count <4000) %>% 
  ggplot(aes(x = complaint_type, y = count, fill = complaint_type)) +
  geom_col() +
  facet_wrap(~ suburb)

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

2. a) Are there any major differences between regions?

2. How do the animal outcomes differ over time?

3. Are there any noticeable anomalies in the data?

It’s impossible to say if there’s different peaks of times of yer per region, as we only have the data for the year as a whole. Unless we look at Suburb rather than region, and break this down by time. However, I think it’d be good to use the general Autralia (nationwide) data.

Proportions overall by region and by outcome. We can see that in general, regions have pretty similar proportions across all outcomes. This won’t be the case for every year though and trends over time, that’s the more interesting stuff.

This is the same thing as before, but for animal type. Same as before, not much discernible difference between the regions.

Next is just the sum ocurrences. We can see the regions that have way more animals injured/taken in by the RSPCA. This could be a population thing, maybe these regions just have more wildlife populations etc., or even a bigger RSPCA presence leading to more cases.

So, this is the data across time (from 1999-2018), with time on the x-axis and number of occurences on the y. Faceted into graphs based on outcome, with the coloured lines representing different regions.

What can we say about this? Well, QLD has had their transfers go way up since about 2009, and in general looks like an increase in cases for them. NSW had a big spike of Currently In Care around 2004. Euthanisation in general seems to be falling out of favour (which is actually true, the RSPCA tries to do it less for various reasons, one actually being the mental health of the staff and another being lots of euthanisations leads to a negative perception of animal shelters)

Exact same as before, but faceted on animal type rather than outcome.

QLD has had a lot more wildlife cases, steadily climbing since about 2006. NSW has had a big drop in cats recently. Dogs too. Horses have gone up every region it looks like. NT has had a huge spike in cats from 2015-2016.

This one has the same x and y, but is faceted based on region. The coloured lines now represent the outcomes.

WA in 2006 had a big increase in euthanizations. VIC has had theirs steadily decreasing over time, this is reflected also in NSW, TAS and ACT.

Same as before with faceting based on region, but now the coloured lines represent the type of animal the RSPCA is dealing with.

NT has had a big spike in dogs and cats in recent years. WA had a big spike in wildlife in 2005, corresponding with their spike in euthanizations - common sense and domain knowledge tell us these are likely highly correlated. VIC’s numbers in general seem to be decreasing over time.

Looking now specifically at euthanizations, faceted by region with the lines representing euthanizations.

Now we no longer facet by region.

OK so cats and dogs have way less euthanisations. Wildlife has gone way up! Horses are too small a sample to be meaningful. Other animals have spiked then decreased.

Looking at rehomes to euthanisations, we see cats have an upward trend with less euthanisations and more rehomes, with dogs, both have gone down. Perhaps with dogs there’s a confounding factor. In fact, it seems like the RSPCA is dealing with less dogs in general.

Euthanizations for wildlife have overall gone up. Queensland is skewing the data, they’ve euthanized way more! Is this necessarily a bad thing? It could indicate a bigger RSPCA presence in the area, and more vigilance in regards to controlling the wildlife.

General thoughts:

big change in 2006 for WA; cyclone Emma March 2006 could’ve released many wild animals leading to high euthanizations (this for 2.c) different territories have different approaches;

4. Business Intelligence – using the insights you have found, can you predict how this might look for the upcoming year?

animal_outcomes %>% 
  drop_na() %>% 
  group_by(region, outcome, year) %>%
  summarize(total_injuries = sum(number_of_occurences)) %>%
  ungroup() %>% 
  filter(outcome == "Euthanized") %>% 
  ggplot(aes(x = year, y = total_injuries, color = outcome)) +
  geom_line() +
  facet_wrap(~region, scales = "free_y") +
  labs(x = "Date", y = "Number of injuries", color = "Outcome")

`summarise()` has grouped output by 'region', 'outcome'. You can override using the `.groups` argument.

# Load required packages
library(forecast)

# Convert the data to a time series object
ts_data <- ts(euthanized_data$total_occurrences, start = c(1999, 1), frequency = 1)

# Fit an ARIMA model to the time series data
arima_model <- auto.arima(ts_data)

# Make forecasts for the next 5 years
forecast_values <- forecast(arima_model, h = 5)

# Plot the forecasted values and the 95% confidence intervals
plot(forecast_values, xlab = "Year", ylab = "Total Euthanized Occurrences", main = "ARIMA Forecast", xlim = c(1999, 2023))

Creating a model for Townsville Total cases:

library(ggthemes)

animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  ggplot() +
  aes(x = date_received, y = count ) +
  geom_line(col = "blue") +
  scale_x_date(date_breaks = "6 months", date_labels = "%b %y") + 
  labs(title = "Townsville Complaint Calls",
       x = "Time",
       y = "Number of Calls") +
  theme_economist() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, family = "mono", face = "bold", size = 7),
        axis.title.y = element_text(vjust = 3.2),
        title = element_text(family = "mono", face = "bold"),
        plot.title = element_text(hjust = 0.4)) +
        annotate("rect", xmin = as.Date("2014-06-01"), xmax = as.Date("2014-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2015-06-01"), xmax = as.Date("2015-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2016-06-01"), xmax = as.Date("2016-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2017-06-01"), xmax = as.Date("2017-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2018-06-01"), xmax = as.Date("2018-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2019-06-01"), xmax = as.Date("2019-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3)

this sorta works:

library(tidyverse)
library(forecast)
library(ggthemes)


# Create a time series from the data
animal_complaints_ts <- animal_complaints %>%
  group_by(date_received) %>%
  summarise(count = n()) %>%
  ts(frequency = 365)

# Forecasting for next 3 years
animal_complaints_fc <- forecast(animal_complaints_ts, h = 3*365)

Error in UseMethod("forecast") : 
  no applicable method for 'forecast' applied to an object of class "c('mts', 'ts', 'matrix')"

library(forecast)
library(ggthemes)

# Convert the date_received column to a Date type

# Aggregate the data to get the count of complaints per day
animal_complaints <- animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  ungroup()

# Convert the data to a time series object
ts_data <- ts(animal_complaints$count, frequency = 365)

# Decompose the time series using STL
ts_stl <- stl(ts_data, s.window = "periodic")

Error in stl(ts_data, s.window = "periodic") : 
  series is not periodic or has less than two periods

# Create a time series with a frequency of 12 (monthly data)
ts_data <- ts(animal_complaints$count, frequency = 12)

# Fit a seasonal ARIMA model to the time series
fit <- auto.arima(ts_data)

# Forecast the next 36 periods (3 years)
forecast_values <- forecast(fit, h = 36)

Error in switch(names(interval), day = "days", hour = "hours", minute = "mins",  : 
  EXPR must be a length 1 vector

ts_data <- ts(data, start = c(2014, 1), frequency = 12)

Error in attr(data, "tsp") <- c(start, end, frequency) : 
  object is not a matrix

library(forecast)

ts_data <- animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  as.data.frame() %>% 
  as.matrix.data.frame()

ts_data <- ts(ts_data, start = c(2014,1), frequency = 12)

fit <- auto.arima(ts_data)

Error in auto.arima(ts_data) : 
  auto.arima can only handle univariate time series

forecast_values <- forecast(fit, h = 36)

Error in switch(names(interval), day = "days", hour = "hours", minute = "mins",  : 
  EXPR must be a length 1 vector